#!/bin/tcsh -ef
####
#
#############################################################################
#                                                                           #
# Title: Diffmap II                                                         #
#                                                                           #
# (C) 2dx.org, GNU Plublic License.                                         #
#                                                                           #
# Created..........: 01/03/2007                                             #
# Last Modification: 01/03/2007                                             #
# Author...........: 2dx.org                                                #
#                                                                           #
#############################################################################
#
# SORTORDER: 71 
#
#
# MANUAL: This is still under the development 
#
#=============================================================================
# Other parameters follow 
#=============================================================================
#
# DISPLAY: MergeResolution
# DISPLAY: zstarrange
# DISPLAY: Calc_from_zstarrange
# DISPLAY: zstarrange_real
# DISPLAY: RESMIN
# DISPLAY: RESMAX
# DISPLAY: merge_res_limit
# DISPLAY: MergeIQMAX
# DISPLAY: MergeHKMAX
# DISPLAY: tempkeep
# DISPLAY: merge_data_type
# DISPLAY: realcell
# DISPLAY: realang
# DISPLAY: ALAT
# DISPLAY: SYM
# DISPLAY: max_amp_correction
# DISPLAY: diff_map1_name
# DISPLAY: diff_map2_name
# DISPLAY: diffmap_cntrs_mode
# DISPLAY: diffmap_cntrs_max
# DISPLAY: diffmap_cntrs_min
# DISPLAY: diffmap_cntrs_step
# DISPLAY: npo_cntrs_step 
# DISPLAY: npo_cntrs_below 
# DISPLAY: diff_map2_dir
# DISPLAY: diffmap_det_significance 
# DISPLAY: diffmap_sig_algorithm 
# DISPLAY: diff_map_confidence
# DISPLAY: diffmap_colormap
# DISPLAY: diffmap_two_by_two
# DISPLAY: diffmap_fixed_plotrange 
# DISPLAY: diffmap_plotrange
# DISPLAY: diffmap_plot_clean 
# DISPLAY: diffmap_plot_scalebar 
# DISPLAY: diffmap_threshold_method
# DISPLAY: diffmap_threshold_global
# DISPLAY: diffmap_var_factor 
#
#$end_local_vars
#
set bin_2dx = ""
set proc_2dx = ""
#
set merge_modus = ""
set zstarwin  = ""
set zstarrange  = ""
set Calc_from_zstarrange = ""
set zstarrange_real = ""
set RESMIN = ""
set RESMAX = ""
set merge_res_limit = ""
set MergeIQMAX = ""
set MergeHKMAX = ""
set merge_data_type = ""
set ILIST = ""
set tempkeep = ""
set realcell = ""
set realang = ""
set ALAT = ""
set SYM = ""
set merge_ref_num = ""
set avrgamphsNUMBER = ""
set avrgamphsRESOL = ""
set max_amp_correction = ""
set diffmap_cntrs_mode = ""
set diffmap_cntrs_max = ""
set diffmap_cntrs_min = ""
set diffmap_cntrs_step = ""
set npo_cntrs_step  = ""
set npo_cntrs_below = ""
#
set lattice = ""
set diff_map2_dir = ""
set diff_map1_name = ""
set diff_map2_name = ""
set diffmap_det_significance = ""
set diffmap_sig_algorithm = "" 
set diff_map_confidence = ""
set diffmap_colormap = ""
set diffmap_two_by_two = ""
set diffmap_fixed_plotrange = ""
set diffmap_plotrange = ""
set diffmap_plot_clean = ""
set diffmap_plot_scalebar = ""
set diffmap_threshold_method = ""
set diffmap_threshold_global = ""
set diffmap_var_factor = ""
#
#$end_vars
#
set scriptname = 2dx_diffmapII
\rm -f LOGS/${scriptname}.results
#
set ccp4_setup = 'y'
source ${proc_2dx}/initialize
#
#################################################################################
${proc_2dx}/linblock "Verifying some parameters"
#################################################################################
#
# This memorizes the current merge directory under the variable "olddir":
set olddir = $PWD
#
echo "The current working directory is" ${olddir}
#
# tried to check for python
#set python_cmd = `where python` || echo "python not found"
#echo "python_cmd = ${python_cmd}"
#if ( ${python_cmd} == "" ) then
#    ${proc_2dx}/protest "python not found!"
#endif
#
if ( ! -d ${diff_map2_dir} ) then
    ${proc_2dx}/protest "The diffmap2 project directory does not exist: ${diff_map2_dir}"
endif
#
set map1 = "map1"
set map2 = "map2"
#
set diffmap_dir = diffmap
if ( ! -d ${diffmap_dir} ) then
    echo "::The directory diffmap does not exist and therefore will be created."
    mkdir ${diffmap_dir}
else
#TODO: uncomment
    rm -rf ${diffmap_dir}
    mkdir ${diffmap_dir}
endif
# get the unit cell size of the 2nd map 
set realcell_map2 = `grep " realcell[ =]"  ${diff_map2_dir}/merge/2dx_merge.cfg  | grep "set" | awk '{print $4}' |  sed 's/"/ /g'`
echo ":real cell ${map2}: ${realcell_map2}"
#
if ( ${Calc_from_zstarrange} == "y" ) then
  set zstarrange_real = `echo ${zstarrange} ${ALAT} | awk '{ s = 1.0 / ( $1 ) } END { print s }'`
  echo "set zstarrange_real = ${zstarrange_real}" >> LOGS/${scriptname}.results
  ${proc_2dx}/linblock "Calculating vertical resolution as ${zstarrange_real} Angstroems."
else
  set zstarrange = `echo ${zstarrange_real} ${ALAT} | awk '{ s = 1.0 / ( $1 ) } END { print s }'`
  echo "set zstarrange = ${zstarrange}" >> LOGS/${scriptname}.results
  ${proc_2dx}/linblock "Calculating zstarrange as ${zstarrange} (with 0.5 = Nyquist resolution)."
endif
#
set zmin = `echo ${zstarrange} | awk '{s = -$1} END {print s}'`
set zminmax = `echo ${zmin},${zstarrange}`
echo zminmax = ${zminmax}
#
#
echo "<<@progress: 1>>"
#
if ( `echo ${RESMIN} ${RESMAX} | awk '{ if ( $1 < $2 ) { s = 1 } else { s = 0 }} END { print s }'` == 1 ) then
  set oldval = ${RESMIN}
  set RESMIN = ${RESMAX}
  set RESMAX = ${oldval}
  ${proc_2dx}/linblock "ERROR: exchanging RESMIN and RESMAX, to RESMIN=${RESMIN}, and RESMAX=${RESMAX}"
  echo "set RESMIN = ${RESMIN}" >> LOGS/${scriptname}.results
  echo "set RESMAX = ${RESMAX}" >> LOGS/${scriptname}.results
endif
#
echo "<<@progress: 5>>"
#
#############################################################################
${proc_2dx}/linblock "sourcing sym2spsgrp_sub.com"
#############################################################################
#
source ${proc_2dx}/2dx_sym2spcgrp_sub.com
#
echo SYM = ${SYM}
echo spcgrp = ${spcgrp}
echo CCP4_SYM = ${CCP4_SYM}
#
############################################################################# 
${proc_2dx}/lin "2dx_merge_makedirs - to create all required subdirectories"
#############################################################################
#
set IS_2DX = yes
source ${proc_2dx}/2dx_merge_makedirs
#
echo "<<@progress: 10>"
#
############################################################################# 
${proc_2dx}/lin "comparing the unit cells size of both maps"
#############################################################################
#
set realcell_a = `echo ${realcell} | sed 's/,/ /g' | awk '{ s = 2*$1  } END { print s }'`
set realcell_b = `echo ${realcell} | sed 's/,/ /g' | awk '{ s = 2*$2  } END { print s }'`
set realcell_map2_a = `echo ${realcell_map2} | sed 's/,/ /g' | awk '{ s = 2*$1  } END { print s }'`
set realcell_map2_b = `echo ${realcell_map2} | sed 's/,/ /g' | awk '{ s = 2*$2  } END { print s }'`
set grid_a = `echo ${realcell} | sed 's/,/ /g' | awk '{ s = 2 * int($1 - ($1 % 4))  } END { print s }'`
set grid_b = `echo ${realcell} | sed 's/,/ /g' | awk '{ s = 2 * int($2 -  ($2 % 4)) } END { print s }'`
set gridsize = "${grid_a} ${grid_b} 20" 
#set realcell_a = `echo ${realcell} | sed 's/,/ /g' | awk '{ s = $1  } END { print s }'`
#set realcell_b = `echo ${realcell} | sed 's/,/ /g' | awk '{ s = $2  } END { print s }'`
#set realcell_map2_a = `echo ${realcell_map2} | sed 's/,/ /g' | awk '{ s = $1  } END { print s }'`
#set realcell_map2_b = `echo ${realcell_map2} | sed 's/,/ /g' | awk '{ s = $2  } END { print s }'`
set unitcell_equal = `echo ${realcell_a} ${realcell_map2_a} ${realcell_b}  ${realcell_map2_b} | awk '{ if( int($1)==int($2) && int($3)==int($4) ) {print 1} else {print 0}}'`
if( ${unitcell_equal} == 0 ) then
    echo ":: The unit cell of both maps are not equal and will be cropped: "
    echo ":: The unit cell map1: ${realcell} map2: ${realcell_map2}"
    if("${diffmap_sig_algorithm}" == "t-Test" ) then
        echo ":: The t-Test can only be used when the unit cell of both maps have the same size"
        echo ":: Therefore the mixed merge method is used"
        set diffmap_sig_algorithm = "Mixed Merge"
    endif
endif
#
if( "${diffmap_sig_algorithm}" == "Welch's t-test" ||  "${diffmap_sig_algorithm}" == "Student's t-test" ) then
    set diffmap_sig_class = "statistical"
else
    set diffmap_sig_class = "empirical"
endif
#
if( "${diffmap_sig_algorithm}" == "Welch's t-test" ) then
    set diffmap_domain = "fourier"
else
    set diffmap_domain = "real"
endif
#
#############################################################################
${proc_2dx}/linblock "2dx_scale_amplitudes.py - scaling amplitudes"
#############################################################################
#TODO: this is also done in 2dx_diffmap scripts, but for the contour to be right we have to do it here as well
set merge_aph1 = "APH/merge.aph"
set merge_aph2 = "${diff_map2_dir}/merge/APH/merge.aph" 
if ( ! -f ${merge_aph1} ) then
    ${proc_2dx}/protest "The merged aph file ${merge_aph1} does not exist!"
endif
if ( ! -f ${merge_aph2} ) then
    ${proc_2dx}/protest "The merged aph file ${merge_aph2} does not exist!"
endif
${app_python} ${bin_2dx}/diffmap/2dx_scale_amplitudes.py ${merge_aph1} 
${app_python} ${bin_2dx}/diffmap/2dx_scale_amplitudes.py ${merge_aph2} 
#############################################################################
${proc_2dx}/linblock "2dx_avrgamphs - only to calculate FOM"
#############################################################################
#
echo "current zminmax = ${zminmax}"
set zmaxlocal = `echo ${ALAT} | awk '{ s = ( 1 / ( 2 * $1 ) ) } END { print s }'`
set zminmaxlocal = `echo -${zmaxlocal},${zmaxlocal}`
# avramphs only works for 2D projection data.
${proc_2dx}/linblock "Symmetry statistics here are only good in 2D."
${proc_2dx}/linblock "Using therefore zminmax=${zminmaxlocal}."
#
# map1
#
\rm -f fort.1
\cp -f ${merge_aph1} fort.1
echo "# IMAGE-IMPORTANT: APH/merge.aph <merge.aph [H,K,Z,A,P,#,IQ,W,Bk,CTF]>" >> LOGS/${scriptname}.results
\rm -f fort.2
\rm -f APH/avrg2D_${map1}.hkl
\rm -f fort.3
\rm -f fort.4
\rm -f TMP444888.tmp
\rm -f TMP444789.tmp
\rm -f LOGS/avramphs_${map1}.table.txt
#
#TODO: check if dir exist clean up
\cp -f APH/merge.aph APH/merge1.aph 
 
${app_python} ${bin_2dx}/diffmap/2dx_diffmap_real.py APH/merge.aph  ${diff_map2_dir}/merge/APH/merge.aph ${zmaxlocal} ${MergeIQMAX} ${max_amp_correction} ${SYM} ${diff_map_confidence}
${bin_2dx}/2dx_avrgamphs.exe << eot > LOGS/2dx_avrgamphs2D_{map1}.log
T
1001,${zminmaxlocal}
${MergeIQMAX}
${avrgamphsNUMBER}
${avrgamphsRESOL}
${realcell} ${realang}
${max_amp_correction}
eot
#
\rm -f fort.1
\mv -f fort.2 APH/avrg2D_${map1}.hkl
\rm -f fort.3
\rm -f fort.4
echo "# IMAGE: LOGS/2dx_avrgamphs2D_{map1}.log <LOG: 2dx_avrgamphs output for 2D run>" >> LOGS/${scriptname}.results
echo "# IMAGE: APH/avrg2D_${map1}.hkl <map1  averaged amp&phs for FOM [H,K,F,P,IQ,FOM]>" >> LOGS/${scriptname}.results
\rm -f TMP444789.tmp 
\rm -f TMP444888.tmp
#
## map2
##
#\rm -f fort.1
#\cp -f ${merge_aph2} fort.1
#echo "# IMAGE-IMPORTANT: APH/merge.aph <merge.aph [H,K,Z,A,P,#,IQ,W,Bk,CTF]>" >> LOGS/${scriptname}.results
#\rm -f fort.2
#\rm -f APH/avrg2D_${map2}.hkl
#\rm -f fort.3
#\rm -f fort.4
#\rm -f TMP444888.tmp
#\rm -f TMP444789.tmp
#\rm -f LOGS/avramphs_${map2}.table.txt
##
##TODO: check if dir exist clean up
#\cp -f ${diff_map2_dir}/merge/APH/merge.aph APH/merge2.aph 
##
#${bin_2dx}/2dx_avrgamphs.exe << eot > LOGS/2dx_avrgamphs2D_{map2}.log
#T
#1001,${zminmaxlocal}
#${MergeIQMAX}
#${avrgamphsNUMBER}
#${avrgamphsRESOL}
#${realcell} ${realang}
#${max_amp_correction}
#eot
##
#\rm -f fort.1
#\mv -f fort.2 APH/avrg2D_${map2}.hkl
#\rm -f fort.3
#\rm -f fort.4
#echo "# IMAGE: LOGS/2dx_avrgamphs2D_{map1}.log <LOG: 2dx_avrgamphs output for 2D run>" >> LOGS/${scriptname}.results
#echo "# IMAGE: APH/avrg2D_${map2}.hkl <map2  averaged amp&phs for FOM [H,K,F,P,IQ,FOM]>" >> LOGS/${scriptname}.results
#\rm -f TMP444789.tmp 
#\rm -f TMP444888.tmp
##
##############################################################################
${proc_2dx}/linblock "2dx_diffmap_sub - avrgamphs with t-test"
#############################################################################
# remove files from previous run
\rm -f APH/avrg2D_sig1.hkl 
\rm -f APH/avrg2D_sig2.hkl 
#
#
echo "<<@progress: 24>>"
#
#TODO: set ${diffmap_nsplits} in GUI
set diffmap_nsplits = 3 
#
if ( "${diffmap_det_significance}" == "y" ) then
    if( "${diffmap_domain}" == "fourier"  && ${unitcell_equal} == 1 ) then
        #Welch's t-test
        rm -rf APH/avrg2D_{map1}.hkl 
        ${app_python} ${bin_2dx}/diffmap/2dx_diffmap_fourier.py APH/merge.aph  ${diff_map2_dir}/merge/APH/merge.aph ${zmaxlocal} ${MergeIQMAX} ${max_amp_correction} ${SYM} ${diff_map_confidence}
    else #real space
        if ( "${diffmap_sig_class}" == "empirical" ) then
            set diffmap_selections_dir = "${diffmap_dir}/selections"
            rm -f ${diffmap_selections_dir}/*.dat
            set diffmap_selection_filebase = "diffmap_split_selection" 
            if( "${diffmap_sig_algorithm}" == "Mixed Merge" ) then
                ${app_python} ${bin_2dx}/diffmap/2dx_merge_split_selections.py --mode mix -n ${diffmap_nsplits} -o ${diffmap_selection_filebase} 2dx_merge_dirfile.dat ${diff_map2_dir}/merge/2dx_merge_dirfile.dat 
            else
                ${app_python} ${bin_2dx}/diffmap/2dx_merge_split_selections.py --mode classic -n ${diffmap_nsplits} -o ${diffmap_selection_filebase} 2dx_merge_dirfile.dat ${diff_map2_dir}/merge/2dx_merge_dirfile.dat 
            endif
            set diffmap_selections = `find ${diffmap_selections_dir} -name ${diffmap_selection_filebase}*.dat`
            echo "${diffmap_selections}"
            set diffmap_merge_output_dir = ${diffmap_selections_dir} 
            foreach selection (${diffmap_selections})
                set diffmap_split_selection = ${selection} 
                set diffmap_selection_stripped = `echo "${diffmap_split_selection}" | awk -F/ '{ print $NF }' `
                set diffmap_merge_output = `echo "${diffmap_selection_stripped}" | cut -d'.' -f1`
                echo ${diffmap_merge_output}
                echo ${diffmap_merge_output_dir}
                source ${proc_2dx}/2dx_diffmap_merge_selection.com
            end
        else # Student's t-Test
            set diffmap_genmaps_dir = "${diffmap_dir}/maps"
            set varmap = "${diffmap_dir}/varmap.mrc"
            rm -rf ${diffmap_genmaps_dir}
            mkdir -p ${diffmap_genmaps_dir}
            set olddir = $PWD
            cd ..
            set basedir = $PWD
            cd ${olddir}
            set diffmap_maps_dir1 = "${diffmap_genmaps_dir}/conformation1" 
            set diffmap_t_test_maps_dir = ${diffmap_maps_dir1} 
            set diffmap_t_test_selection = "2dx_merge_dirfile.dat" 
            source  ${bin_2dx}/diffmap/2dx_diffmap_t-test.com
            set diffmap_maps_dir2 = "${diffmap_genmaps_dir}/conformation2" 
            set diffmap_t_test_maps_dir = ${diffmap_maps_dir2} 
            set basedir = ${diff_map2_dir}
            set diffmap_t_test_selection = "${diff_map2_dir}/merge/2dx_merge_dirfile.dat"
            source  ${bin_2dx}/diffmap/2dx_diffmap_t-test.com
            ${app_python} ${bin_2dx}/diffmap/2dx_diffmap_t-test.py --confidence ${diff_map_confidence} -o ${varmap} ${diffmap_maps_dir1} ${diffmap_maps_dir2}
            #redo merge
            ${proc_2dx}/linblock "remerging map 1"
            set diffmap_split_selection = "2dx_merge_dirfile.dat"
            set diffmap_merge_output = merge 
            set diffmap_merge_output_dir = APH 
            source ${proc_2dx}/2dx_diffmap_merge_selection.com
            # redo merge for map2
            cd ${diff_map2_dir}/merge
            echo ":: current working dir $PWD"
            ${proc_2dx}/linblock "remerging map 2"
            set diffmap_split_selection = "2dx_merge_dirfile.dat"
            set diffmap_merge_output = merge 
            set diffmap_merge_output_dir = "${diff_map2_dir}/merge/APH" 
            source ${proc_2dx}/2dx_diffmap_merge_selection.com
            cd ${olddir}
        endif 
        ${app_python} ${bin_2dx}/diffmap/2dx_diffmap_real.py  APH/merge.aph  ${diff_map2_dir}/merge/APH/merge.aph ${zmaxlocal} ${MergeIQMAX} ${max_amp_correction} ${SYM} ${diff_map_confidence}
    endif
else
    ${app_python} ${bin_2dx}/diffmap/2dx_diffmap_real.py APH/merge.aph  ${diff_map2_dir}/merge/APH/merge.aph ${zmaxlocal} ${MergeIQMAX} ${max_amp_correction} ${SYM} ${diff_map_confidence}
    #TODO: change all the image names for the non significant case
endif
#
echo "# IMAGE-IMPORTANT: APH/avrg2D_sig1.hkl <map1 significant reflections [H,K,Z,A,P,FOM]>" >> LOGS/${scriptname}.results
echo "# IMAGE-IMPORTANT: APH/avrg2D_sig2.hkl <map2 significant reflections [H,K,Z,A,P,FOM]>" >> LOGS/${scriptname}.results
#
#############################################################################
${proc_2dx}/linblock "2dx_centric2 - to correct phases to 0 or 180 for 2D run"
#############################################################################  
#
#TODO: only if domain is fourier 
# map1
#
\rm -f APH/centric2D_${map1}.hkl
\rm -f APH/centric2D_${map1}.hk
#
${bin_2dx}/2dx_centric2.exe << eot
APH/avrg2D_${map1}.hkl
APH/centric2D_${map1}.hkl
APH/centric2D_${map1}.hk
${realcell},${realang}
${RESMIN},${RESMAX}
${SYM}
eot
#
if ( ! -e APH/centric2D_${map1}.hkl ) then
  ${proc_2dx}/protest "ERROR occured."
endif
# sigmap1
#
\rm -f APH/centric2D_${map1}_sig.hkl
\rm -f APH/centric2D_${map1}_sig.hk
#
if ( ! -e APH/avrg2D_sig1.hkl ) then
  ${proc_2dx}/protest "avrg2D_sig1.hkl does not exist."
else
  echo "avrg2D_sig1.hkl exists."
endif
#
${bin_2dx}/2dx_centric2.exe << eot
APH/avrg2D_sig1.hkl
APH/centric2D_${map1}_sig.hkl
APH/centric2D_${map1}_sig.hk
${realcell},${realang}
${RESMIN},${RESMAX}
${SYM}
eot
#
if ( ! -e APH/centric2D_${map1}_sig.hkl ) then
  ${proc_2dx}/protest "ERROR occured."
endif
#
# sigmap2
#
\rm -f APH/centric2D_${map2}_sig.hkl
\rm -f APH/centric2D_${map2}_sig.hk
#
${bin_2dx}/2dx_centric2.exe << eot
APH/avrg2D_sig2.hkl
APH/centric2D_${map2}_sig.hkl
APH/centric2D_${map2}_sig.hk
${realcell_map2},${realang}
${RESMIN},${RESMAX}
${SYM}
eot
#
echo "<<@progress: 30>>"
#
if ( ! -e APH/centric2D_${map2}_sig.hkl ) then
  ${proc_2dx}/protest "ERROR occured."
endif
#
# This is not used:
\rm -f APH/centric2D_{map1}.hk
\rm -f APH/centric2D_{map2}.hk
#
#if ( ${tempkeep} == "y" ) then
  echo "# IMAGE: APH/centric2D_${map1}_sig.hkl <map1 after CENTRIC for 2D [H,K,L,F,P,FOM]>" >> LOGS/${scriptname}.results
  echo "# IMAGE: APH/centric2D_${map2}_sig.hkl <map2 after CENTRIC for 2D [H,K,L,F,P,FOM]>" >> LOGS/${scriptname}.results
#endif
#
echo "<<@progress: 35>>"
#
#############################################################################
${proc_2dx}/linblock "2dx_hklsym4 - to apply symmetry to APH file for 2D run"
#############################################################################  
# TODO: Replace this and symmetrize already earlier
#
\rm -f APH/sym2D_{map1}.hkl
\rm -f APH/sym_nosort2D_{map1}.hkl
\rm -f APH/sym_sort2D_{map1}.hkl
\rm -f APH/sym_noheader2D_{map1}.hkl
#
# Set isig to 3, for NO SIGF BUT SET SIGF to 1.0
set isig = 3
#
# map1
${bin_2dx}/2dx_hklsym4.exe << eot
APH/centric2D_${map1}.hkl
APH/sym_nosort2D_${map1}.hkl
APH/sym2D_${map1}.hkl
${spcgrp}
1
${isig}
0     ! write out full p1 plane
0     ! do not write out negative L values
eot
#
# LABOUT H K L F PHI FOM SIGF
# CTYPOUT H H H F P W Q
#
#echo "# IMAGE: APH/sym2D_${map1}.hkl <map1 APH after symmetrization [H,K,L,F,P,FOM,1.0]>" >> LOGS/${scriptname}.results
#
if ( ! -e APH/sym2D_${map1}.hkl ) then
  ${proc_2dx}/protest "ERROR occured."
endif
#
# sigmap1
\rm -f APH/sym2D_{map1}_sig.hkl
\rm -f APH/sym_nosort2D_{map1}_sig.hkl
\rm -f APH/sym_sort2D_{map1}_sig.hkl
\rm -f APH/sym_noheader2D_{map1}_sig.hkl
#
${bin_2dx}/2dx_hklsym4.exe << eot
APH/centric2D_${map1}_sig.hkl
APH/sym_nosort2D_${map1}_sig.hkl
APH/sym2D_${map1}_sig.hkl
${spcgrp}
1
${isig}
0     ! write out full p1 plane
0     ! do not write out negative L values
eot
#
echo "# IMAGE: APH/sym2D_${map1}_sig.hkl <map1 APH after symmetrization [H,K,L,F,P,FOM,1.0]>" >> LOGS/${scriptname}.results
#
if ( ! -e APH/sym2D_${map1}_sig.hkl ) then
  ${proc_2dx}/protest "ERROR occured."
endif
#
# sigmap2
\rm -f APH/sym2D_{map2}_sig.hkl
\rm -f APH/sym_nosort2D_{map2}_sig.hkl
\rm -f APH/sym_sort2D_{map2}_sig.hkl
\rm -f APH/sym_noheader2D_{map2}_sig.hkl
#
${bin_2dx}/2dx_hklsym4.exe << eot
APH/centric2D_${map2}_sig.hkl
APH/sym_nosort2D_${map2}_sig.hkl
APH/sym2D_${map2}_sig.hkl
${spcgrp}
1
${isig}
0     ! write out full p1 plane
0     ! do not write out negative L values
eot
#
#
echo "# IMAGE: APH/sym2D_${map2}_sig.hkl <map2 APH after symmetrization [H,K,L,F,P,FOM,1.0]>" >> LOGS/${scriptname}.results
#
if ( ! -e APH/sym2D_${map2}_sig.hkl ) then
  ${proc_2dx}/protest "ERROR occured."
endif
#
############################################################################
${proc_2dx}/linblock "f2mtz - to transform APH file into MTZ file for 2D run"
#############################################################################  
#
set LABOUTval = "H K L F PHI FOM"
set CTYPOUTval = "H H H F P W"
#
# map1
set infile = APH/sym2D_${map1}.hkl
set outfile = SCRATCH/${map1}_MRClefthanded.mtz
\rm -f SCRATCH/${map1}_MRClefthanded.mtz
#
${bin_ccp4}/f2mtz hklin ${infile} hklout ${outfile} << eof
TITLE  Map, Symmetry=${CCP4_SYM}, ${map1} , ${date}
#CELL ${realcell_a} ${realcell_b} ${ALAT} 90.0 90.0 ${realang}
CELL ${realcell} ${ALAT} 90.0 90.0 ${realang}
SYMMETRY ${CCP4_SYM} 
LABOUT ${LABOUTval}
CTYPOUT ${CTYPOUTval}
FILE ${infile}
SKIP 0
END
eof
#
\rm -f ${map1}_MRClefthanded.mtz 
${bin_ccp4}/sftools << eot
read ${outfile} 
merge
expand
write ${map1}_MRClefthanded.mtz
end
eot
#TODO: This is still just a test
#${bin_ccp4}/sftools << eot
#read ${map1}_MRClefthanded.mtz
#calc P col PHI = col 2 3.14159 +
#delete col 2
#write ${map1}_shifted_MRClefthanded.mtz
#END
#eot
#
#\rm -f ${map1}_MRClefthanded.mtz 
#cp ${map1}_shifted_MRClefthanded.mtz ${map1}_MRClefthanded.mtz
# sigmap1
set infile = APH/sym2D_${map1}_sig.hkl
set outfile = SCRATCH/${map1}_sig_MRClefthanded.mtz
\rm -f SCRATCH/${map1}_sig_MRClefthanded.mtz
#
${bin_ccp4}/f2mtz hklin ${infile} hklout ${outfile} << eof
TITLE  Map, Symmetry=${CCP4_SYM}, ${diff_map1_name} , ${date}
#CELL ${realcell_a} ${realcell_b} ${ALAT} 90.0 90.0 ${realang}
CELL ${realcell} ${ALAT} 90.0 90.0 ${realang}
SYMMETRY ${CCP4_SYM} 
LABOUT ${LABOUTval}
CTYPOUT ${CTYPOUTval}
FILE ${infile}
SKIP 0
END
eof
#
\rm -f ${map1}_sig_MRClefthanded.mtz 
${bin_ccp4}/sftools << eot
read ${outfile} 
merge
expand
write ${map1}_sig_MRClefthanded.mtz
end
eot
#
# sigmap2
set infile = APH/sym2D_${map2}_sig.hkl
set outfile = SCRATCH/${map2}_sig_MRClefthanded.mtz
\rm -f SCRATCH/${map2}_sig_MRClefthanded.mtz
#
${bin_ccp4}/f2mtz hklin ${infile} hklout ${outfile} << eof
TITLE  Map, Symmetry=${CCP4_SYM}, ${diff_map2_name} , ${date}
#CELL ${realcell_map2_a} ${realcell_map2_b} ${ALAT} 90.0 90.0 ${realang}
CELL ${realcell_map2} ${ALAT} 90.0 90.0 ${realang}
SYMMETRY ${CCP4_SYM}
LABOUT ${LABOUTval}
CTYPOUT ${CTYPOUTval}
FILE ${infile}
SKIP 0
END
eof
#
\rm -f ${map2}_sig_MRClefthanded.mtz 
${bin_ccp4}/sftools << eot
read ${outfile} 
merge
expand
write ${map2}_sig_MRClefthanded.mtz
end
eot
#
echo "# IMAGE: ${map1}_sig_MRClefthanded.mtz <MTZ: map 1 merged full 2D data>" >> LOGS/${scriptname}.results
echo "# IMAGE: ${map2}_sig_MRClefthanded.mtz <MTZ: map 2 merged full 2D data>" >> LOGS/${scriptname}.results
#
##############################################################################
${proc_2dx}/linblock  "Inverse fft to get map" 
##############################################################################
#TODO: have it determined by the unit cell size
#set grid_a = `echo ${realcell_a} | sed 's/,/ /g' | awk '{ s = 10 * int($1)  } END { print s }'`
#set grid_b = `echo ${realcell_b} | sed 's/,/ /g' | awk '{ s = 10 * int($1) } END { print s }'`
set xlim1 = `echo ${grid_a} | awk '{ s = $1 - 1} END { print s }'`
set ylim1 = `echo ${grid_b} | awk '{ s = $1 - 1} END { print s }'`
echo  ":: The unit cell was cropped to ${grid_a}, ${grid_b}, with the limit ${xlim1} and ${ylim1}" 

#set grid2_a = `echo ${realcell_map2_a} | sed 's/,/ /g' | awk '{ s = 10 * int($1)  } END { print s }'`
#set grid2_b = `echo ${realcell_map2_b} | sed 's/,/ /g' | awk '{ s = 10 * int($1)  } END { print s }'`
set grid2_a = `echo ${realcell_map2} | sed 's/,/ /g' | awk '{ s = 2 * int($1 -  ($1 % 4))  } END { print s }'`
set grid2_b = `echo ${realcell_map2} | sed 's/,/ /g' | awk '{ s = 2 * int($2 -  ($2 % 4))  } END { print s }'`
set xlim2 = `echo ${grid2_a} | awk '{ s = $1 - 1} END { print s }'`
set ylim2 = `echo ${grid2_b} | awk '{ s = $1 - 1} END { print s }'`
#echo  ": The unit cell was cropped to ${grid2_a}, ${grid2_b}, with the limit ${xlim2} and ${ylim2}" 

# THIS WAS ONLY USED FOR CCP4 CROPING
# which is not used anymore
#set realcell_max_a = `echo ${realcell_a} ${realcell_map2_a} | awk '{ if( $1 > $2 ) {print $1} else {print $2} }'`
#set realcell_max_b = `echo ${realcell_b} ${realcell_map2_b} | awk '{ if( $1 > $2 ) {print $1} else {print $2} }'`
set realcell_min_a = `echo ${realcell_a} ${realcell_map2_a} | awk '{ if( $1 < $2 ) {print $1} else {print $2} }'`
set realcell_min_b = `echo ${realcell_b} ${realcell_map2_b} | awk '{ if( $1 < $2 ) {print $1} else {print $2} }'`
#set realcell_max = "${realcell_max_a} ${realcell_max_b}"
set grid_min_a = `echo ${realcell_min_a} | awk '{ s = int($1)+1 } END { print s }'`
set grid_min_b = `echo ${realcell_min_b} | awk '{ s = int($1)+1 } END { print s }'`
#set xlim_min = `echo ${grid_min_a} | awk '{ s = $1 - 1} END { print s }'`
#set ylim_min = `echo ${grid_min_b} | awk '{ s = $1 - 1} END { print s }'`
set xlim_2x = `echo ${grid_min_a} | awk '{ s = 2*$1 - 1} END { print s }'`
set ylim_2x = `echo ${grid_min_b} | awk '{ s = 2*$1 - 1} END { print s }'`
#
# map1
\rm -f ${map1}.map
#
${bin_ccp4}/fft HKLIN  ${map1}_MRClefthanded.mtz MAPOUT  ${map1}.map \
    << eot
LABIN F1=F PHI=PHI W=FOM
TITLE "${RESMAX} EM map for ${map1}"
PROJECTION
AXIS X Y Z
SCALE F1 1.0
RESOLUTION ${RESMIN} ${RESMAX}
#GRID SAMPLE 20 
GRID  ${grid_a} ${grid_b} 20 
XYZLIM ASU
END
eot
#
# sigmap1
\rm -f ${map1}_sig.map
#
${bin_ccp4}/fft HKLIN  ${map1}_sig_MRClefthanded.mtz MAPOUT  ${map1}_sig.map \
    << eot
LABIN F1=F PHI=PHI W=FOM
TITLE "${RESMAX} EM map for ${diff_map1_name}"
PROJECTION
AXIS X Y Z
SCALE F1 1.0
RESOLUTION ${RESMIN} ${RESMAX}
#GRID SAMPLE 20 
GRID  ${grid_a} ${grid_b} 20 
XYZLIM ASU 
END
eot
#
# sigmap2
\rm -f ${map2}_sig.map
#
${bin_ccp4}/fft HKLIN  ${map2}_sig_MRClefthanded.mtz MAPOUT  ${map2}_sig.map \
    << eot
LABIN F1=F PHI=PHI W=FOM
TITLE "${RESMAX} EM map for ${diff_map2_name}"
PROJECTION
AXIS X Y Z
SCALE F1 1.0
RESOLUTION ${RESMIN} ${RESMAX}
#GRID SAMPLE 20 
GRID  ${grid2_a} ${grid2_b} 20
#XYZLIM 0 ${xlim2} 0 ${ylim2} 0 0  
XYZLIM ASU 
END
eot
#

if( "${diffmap_sig_class}" == "empirical" ) then
    if ( "${diffmap_det_significance}" == "y") then
        set merge2map_output_dir = ${diffmap_selections_dir}
        set merged_selections = `find diffmap -name "*.aph"`
        foreach merge_aph ($merged_selections)
            #############################################################################
            ${proc_2dx}/linblock "Scaling amplitudes of ${merge_aph}"
            #############################################################################
            ${app_python} ${bin_2dx}/diffmap/2dx_scale_amplitudes.py ${merge_aph} 
            set postfix = ` echo ${merge_aph} | awk -F /  '{print $NF}' | sed s/.aph//`
            source ${proc_2dx}/2dx_from_merge2map.com
            set map_postfix = "${merge2map_output_dir}/map_${postfix}"
            if ( ! -e ${map_postfix}_MRClefthanded.mtz ) then
                ${proc_2dx}/protest "${map_postfix}_MRClefthanded.mtz does not exist."
            endif
            \rm -f ${map_postfix}.map 
            #
            ${bin_ccp4}/fft HKLIN ${map_postfix}_MRClefthanded.mtz  MAPOUT ${map_postfix}.map << eot
            LABIN F1=F PHI=PHI W=FOM
            TITLE "${RESMAX} EM map for ${postfix}"
            PROJECTION
            AXIS X Y Z
            SCALE F1 1.0
            RESOLUTION ${RESMIN} ${RESMAX}
            GRID  ${grid_a} ${grid_b} 20 
            XYZLIM ASU
            END
eot
            if( "${diffmap_two_by_two}" == "y" ) then
                ${bin_ccp4}/mapmask mapin ${map_postfix}.map mapout ${map_postfix}.map <<eof
                XYZLIM 0 ${xlim_2x} 0 ${ylim_2x} 0 0
                END
eof
            endif
            #
            \rm -f ${map_postfix}.mrc
            #
${bin_2dx}/labelh.exe << eot
${map_postfix}.map 
2
${map_postfix}.mrc
1,0
0
eot
        end
#       END FOREACH
    endif
endif
#
##############################################################################
${proc_2dx}/linblock  "2dx_variationmap - merge all variation maps to one map" 
##############################################################################
if( "${diffmap_det_significance}" == "y" && "${diffmap_sig_class}" == "empirical" ) then
    set varmap = "${diffmap_dir}/varmap.mrc"
    set variationmap_options =  "-m ${diffmap_threshold_method} -o ${varmap}"
    if( "${diffmap_threshold_global}" == "y" ) then
        set variationmap_options = "${variationmap_options} -g" 
    endif
    ${app_python} ${bin_2dx}/diffmap/2dx_variationmap.py ${diffmap_selections_dir} ${variationmap_options} 
endif
#
#
if( "${diffmap_two_by_two}" == "y" ) then
    ${bin_ccp4}/mapmask mapin ${map1}.map mapout ${map1}.map <<eof
        XYZLIM 0 ${xlim_2x} 0 ${ylim_2x} 0 0
        END
eof

    ${bin_ccp4}/mapmask mapin ${map1}_sig.map mapout ${map1}_sig.map <<eof
        XYZLIM 0 ${xlim_2x} 0 ${ylim_2x} 0 0
        END
eof

    ${bin_ccp4}/mapmask mapin ${map2}_sig.map mapout ${map2}_sig.map <<eof
        XYZLIM 0 ${xlim_2x} 0 ${ylim_2x} 0 0
        END
eof
endif

echo "<<@progress: 30>>"
##############################################################################
${proc_2dx}/linblock  "npo - to create a line plot for both maps" 
##############################################################################
# map1
${bin_ccp4}/npo  MAPIN ${map1}.map PLOT ${map1}.plot <<EOF
TITLE   ${map1}
MAP SCALE 0.5
CONTRS SIG -3.0  to 3.0 BY ${npo_cntrs_step} 
MODE BELOW 0 DASHED 1 0.25 0
MODE BELOW ${npo_cntrs_below}  DASHED 1 0.25 0
XYZLIM 0 ${xlim_2x} 0 ${ylim_2x} 0 0 
SECTS 0 0
PLOT 
EOF
#
# sigmap1
${bin_ccp4}/npo  MAPIN ${map1}_sig.map PLOT ${map1}_sig.plot <<EOF
TITLE   ${diff_map1_name}
MAP SCALE 0.5
CONTRS SIG -3.0  to 3.0 BY ${npo_cntrs_step} 
MODE BELOW 0 DASHED 1 0.25 0
MODE BELOW ${npo_cntrs_below}  DASHED 1 0.25 0
XYZLIM 0 ${xlim_2x} 0 ${ylim_2x} 0 0 
SECTS 0 0
PLOT 
EOF
# 
# sigmap2
${bin_ccp4}/npo  MAPIN ${map2}_sig.map PLOT ${map2}_sig.plot <<EOF
TITLE   ${diff_map2_name}
MAP SCALE 0.5
CONTRS SIG -3.0  to 3.0 BY ${npo_cntrs_step} 
MODE BELOW 0 DASHED 1 0.25 0
MODE BELOW ${npo_cntrs_below}  DASHED 1 0.25 0
XYZLIM 0 ${xlim_2x} 0 ${ylim_2x} 0 0 
SECTS 0 0
PLOT 
EOF
#
#
echo "<<@progress: 35>>"
#############################################################################
${proc_2dx}/linblock "laserplot - to create contour plots"
#############################################################################
#
# map1
\rm -f PS/${map1}.ps
${bin_2dx}/laserplot.exe -outputfile=PS/${map1}.ps ${map1}.plot
\rm -f ${map1}.plot
echo "# IMAGE-IMPORTANT: PS/${map1}.ps <PS: map 1>"  >> LOGS/${scriptname}.results
#
# sigmap1
\rm -f PS/${map1}_sig.ps
${bin_2dx}/laserplot.exe -outputfile=PS/${map1}_sig.ps ${map1}_sig.plot
\rm -f ${map1}_sig.plot
echo "# IMAGE-IMPORTANT: PS/${map1}_sig.ps <PS: significant map 1>"  >> LOGS/${scriptname}.results
#
# sigmap2
\rm -f PS/${map2}_sig.ps
${bin_2dx}/laserplot.exe -outputfile=PS/${map2}_sig.ps ${map2}_sig.plot
\rm -f ${map2}_sig.plot
echo "# IMAGE-IMPORTANT: PS/${map2}_sig.ps <PS: significant map 2>"  >> LOGS/${scriptname}.results
#############################################################################
${proc_2dx}/linblock "LABEL - to create a clean MRC file format instead of CCP4"
#############################################################################
#
# map1
\rm -f ${map1}.mrc
#
${bin_2dx}/labelh.exe << eot
${map1}.map 
2
${map1}.mrc
1,0
0
eot
#
echo "# IMAGE-IMPORTANT: ${map1}.mrc <map 1>" >> LOGS/${scriptname}.results
# sigmap1
\rm -f ${map1}_sig.mrc
#
${bin_2dx}/labelh.exe << eot
${map1}_sig.map 
2
${map1}_sig.mrc
1,0
0
eot
#
echo "# IMAGE-IMPORTANT: ${map1}_sig.mrc <significant map 1>" >> LOGS/${scriptname}.results
#
# sigmap2
\rm -f ${map2}_sig.mrc
#
${bin_2dx}/labelh.exe << eot
${map2}_sig.map 
2
${map2}_sig.mrc 
1,0
0
eot
#
echo "# IMAGE-IMPORTANT: ${map2}_sig.mrc <significant map 2>" >> LOGS/${scriptname}.results
#
#############################################################################
${proc_2dx}/linblock "labelh - to transform the map into BYTE format with automatic scaling"
#############################################################################

${bin_2dx}/labelh.exe << eot
${map1}.mrc
-3
SCRATCH/${map1}.mrc
eot
#
${bin_2dx}/labelh.exe << eot
${map1}_sig.mrc
-3
SCRATCH/${map1}_sig.mrc
eot
#
${bin_2dx}/labelh.exe << eot
${map2}_sig.mrc
-3
SCRATCH/${map2}_sig.mrc
eot


#TODO: let this be a user variable
#set diffmap_shift180 = 'y'
#
# SHIFT by 180 degrees in u direction for plane groups
set shift_180_x = "" 
if ( ${SYM} == 'p4212' ||  ${SYM} == 'p121_b' || ${SYM} == 'p2221b' ) then
    set shift_180_x = "-s"
endif
#
#############################################################################
${proc_2dx}/linblock "2dx_plot_diffmap - to plot the diffmap"
#############################################################################
#
set diffmap_plot_options =   "--map1-name ${diff_map1_name} --map2-name ${diff_map2_name} --colormap ${diffmap_colormap} ${shift_180_x} " 
if ( ${diffmap_det_significance} == "y" ) then
    if ( ${diffmap_domain} == "real" ) then
        set diffmap_plot_options = "--varmap ${varmap} --variation-factor ${diffmap_var_factor} ${diffmap_plot_options}"
    else
        # t-Welch
        set diffmap_plot_options = "--contour-map ${map1}.mrc --map1-name ${diff_map1_name} --map2-name ${diff_map2_name} --colormap ${diffmap_colormap} ${shift_180_x} "
    endif 
endif
if ( ${diffmap_fixed_plotrange} == "y" ) then
    set diffmap_plot_options = "${diffmap_plot_options} --plot-range ${diffmap_plotrange}"
endif
if ( ${diffmap_plot_clean} == "y" ) then
    set diffmap_plot_options = "${diffmap_plot_options} --plot-map-only "
endif
if ( ${diffmap_plot_scalebar} == "y" ) then
    set diffmap_plot_options = "${diffmap_plot_options} --plot-scalebar "
endif
${app_python} ${bin_2dx}/diffmap/2dx_plot_diffmap.py ${diffmap_plot_options}  ${map1}_sig.mrc ${map2}_sig.mrc  
#
# Get the diffmape name
set diffmap_name = `find diffmap -d 1 -name "diffmap*.pdf" | head -1`
if ( ! -z ${diffmap_name} ) then
    echo "# IMAGE-IMPORTANT: ${diffmap_name} <difference map>" >> LOGS/${scriptname}.results
endif
#
#############################################################################
${proc_2dx}/linblock "2dx_merge normal end."
#############################################################################
#
echo "<<@progress: 100>>"
#
